{
"cells": [
{
"cell_type": "markdown",
"id": "d8bb6046",
"metadata": {},
"source": [
"# Optical Observations of Satellites\n",
"\n",
"Ground-based telescopes can observe sunlit satellites against the night sky, measuring their angular position (right ascension and declination). These line-of-sight measurements can then refine the satellite's orbit. This technique is particularly valuable for satellites in MEO and GEO orbits, which are visible for long periods at night and lack easy access to GPS signals. Organizations like the US Space Force and commercial companies such as [ExoAnalytic](https://exoanalytic.com/space-intelligence/) routinely use optical tracking to maintain satellite catalogs.\n",
"\n",
"## Approach\n",
"\n",
"This example uses a **batch least-squares** (Gauss-Newton) technique to update an initial state estimate by minimizing the difference between observed and predicted line-of-sight unit vectors:\n",
"\n",
"$$\n",
"\\hat{\\mathbf x}_0\n",
"=\n",
"\\arg\\min_{\\mathbf x_0\\in\\mathbb R^6}\n",
"\\sum_{k=1}^{N}\n",
"\\left\\|\n",
"E_k\\Big(\\hat{\\mathbf u}_m(\\alpha_k,\\delta_k)-\\hat{\\mathbf u}(\\phi(t_k,t_0;\\mathbf x_0),\\mathbf r_{0}(t_k))\\Big)\n",
"\\right\\|_{R^{-1}}^{2}\n",
"$$\n",
"\n",
"where:\n",
"\n",
"- $\\hat{\\mathbf u}_m(\\alpha,\\delta) = [\\cos\\delta\\cos\\alpha,\\; \\cos\\delta\\sin\\alpha,\\; \\sin\\delta]^T$ is the measured unit vector from right ascension $\\alpha$ and declination $\\delta$\n",
"- $\\hat{\\mathbf u}(\\mathbf x_k, \\mathbf r_0) = \\frac{\\mathbf r_k - \\mathbf r_0}{\\|\\mathbf r_k - \\mathbf r_0\\|}$ is the predicted unit vector from observer at $\\mathbf r_0$ to satellite at $\\mathbf r_k$\n",
"- $\\mathbf x_k = \\phi(t_k, t_0; \\mathbf x_0)$ is the propagated state at time $t_k$ given initial state $\\mathbf x_0$\n",
"- $E_k$ projects the 3D residual into the local tangent plane (RA/Dec directions) to keep the measurement space 2D and avoid singularities\n",
"\n",
"The state transition matrix $\\Phi$ from the high-precision propagator maps initial state perturbations to perturbations at each measurement time, forming the Jacobian of the linearized system."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a53d5d64",
"metadata": {},
"outputs": [],
"source": "import satkit as sk\nimport numpy as np\nimport math as m\n\n# Take our observer to be the Space Surveillance Telescope (SST)\n# in Exmouth, Western Australia\nsst_lat = -22.675\nsst_lon = 114.235\nsst_alt = 0.5 # km\nsst = sk.itrfcoord(latitude_deg=sst_lat, longitude_deg=sst_lon, alt_m=sst_alt * 1e3)\n\n# Fetch the current TLE for INTELSAT 19 (NORAD 38356) directly from CelesTrak.\n# IS-19 is a geostationary communications satellite at 166°E longitude,\n# visible from the SST.\niss19_url = \"https://celestrak.org/NORAD/elements/gp.php?CATNR=38356&FORMAT=TLE\"\ntle = sk.TLE.from_url(iss19_url)\nepoch = tle.epoch\n\n# Lets make the initial \"true\" state the tle-predicted state at epoch\npteme, vteme = sk.sgp4(tle, epoch)\nq = sk.frametransform.qteme2gcrf(epoch)\npgcrf = q * pteme\nvgcrf = q * vteme\nstate0 = np.hstack((pgcrf, vgcrf))\nprint(epoch)\nbegin_time = epoch\nend_time = epoch + sk.duration(days=2)\nsettings = sk.propsettings()\nsettings.precompute_terms(begin_time, end_time)\n# True state over 1 day\ntruth = sk.propagate(state0, begin_time, end_time, propsettings=settings)\n\n# Assume position knowledge to 10 km, velocity to 5 cm/s\npos_noise = 10e3\nvel_noise = 0.05\n\n# initial state estimate\nstate0_est = state0 + np.random.normal(scale=[pos_noise] * 3 + [vel_noise] * 3)\nstate0_est_prior = state0_est.copy()\n\n# Measurements taken once every 5 minutes over a 7-hour window starting\n# 9 hours after the TLE epoch (approximating a single night of observation).\nsample_start = epoch + sk.duration(hours=9)\nsample_end = sample_start + sk.duration(hours=7)\nsample_times = [sample_start + sk.duration(minutes=5) * i for i in range(7 * 60 // 5)]\n\n# Truth at first sample time\nstate0_truth = truth.interp(sample_times[0])\n\n# the optical measurement is a line-of-sight vector from observer to satellite in the inertial frame\nsst_pos = [sk.frametransform.qitrf2gcrf(t) * sst.vector for t in sample_times]\nsat_state = [truth.interp(t) for t in sample_times]\nlos_meas = [sat_state[idx][0:3] - sst_pos[idx] for idx in range(len(sample_times))]\n# Normalize the line-of-sight measurements to get unit vectors\nlos_meas = [los / np.linalg.norm(los) for los in los_meas]\n\n# add noise to the measurements. Assume 30 microradians of angular noise\n# Note: I made up 30 microradians (6 arcseconds).\nang_noise = 30e-6\n\n\ndef add_ang_noise(u, ang_noise):\n # add noise in a random direction perpendicular to the line of sight\n # first get a random vector\n rand_vec = np.random.normal(size=3)\n # make it perpendicular to the line of sight\n rand_vec -= rand_vec.dot(u) * u\n rand_vec /= np.linalg.norm(rand_vec)\n # rotate the line of sight by the noise angle in the direction of the random vector\n theta = np.random.normal(scale=ang_noise)\n u_noisy = (\n u * m.cos(theta)\n + np.cross(rand_vec, u) * m.sin(theta)\n + rand_vec * rand_vec.dot(u) * (1 - m.cos(theta))\n )\n return u_noisy\n\n\nlos_meas_noisy = [add_ang_noise(los, ang_noise) for los in los_meas]\n\n\n# Now, make up an initial state estimate that we are going to refine with measurements\nstate0_est = truth.interp(sample_times[0]) + np.random.normal(\n scale=[pos_noise] * 3 + [vel_noise] * 3\n)\n\n# Our initial state error .. record for later comparison\nstate0_est_prior = state0_est.copy()\n\n\ndef unit_vector_jacobian(rho):\n \"\"\"\n Jacobian of unit vector with respect to the original vector\n \"\"\"\n rho_norm = np.linalg.norm(rho)\n u = rho / rho_norm\n return (np.eye(3) - np.outer(u, u)) / rho_norm\n\n\ndef tangent_basis(u):\n \"\"\"\n RA/Dec tangent basis at predicted direction\n \"\"\"\n ux, uy, uz = u\n alpha = np.arctan2(uy, ux)\n delta = np.arcsin(uz)\n\n e_alpha = np.array([-np.sin(alpha), np.cos(alpha), 0.0])\n e_delta = np.array(\n [-np.sin(delta) * np.cos(alpha), -np.sin(delta) * np.sin(alpha), np.cos(delta)]\n )\n return np.vstack((e_alpha, e_delta))\n\n\n# OK, lets do it!\n# 3 iterations of Gauss-Newton least squares\nfor n in range(0, 3):\n print(f\"Iteration {n}\")\n\n # Propagate state to get state transition matrix at each time\n res = sk.propagate(\n state0_est,\n sample_times[0],\n sample_times[-1],\n output_phi=True,\n propsettings=settings,\n )\n [state, phi] = zip(*[res.interp(t, output_phi=True) for t in sample_times])\n\n residuals = []\n J_list = []\n\n for idx, t in enumerate(sample_times):\n # Get the predicted line of sight measurement at this time\n p_sat = state[idx][0:3]\n p_obs = sst_pos[idx]\n los_vec = p_sat - p_obs\n los_pred = los_vec / np.linalg.norm(los_vec)\n\n # Get the tangent basis in predicted direction\n E = tangent_basis(los_pred)\n\n # Get residuals\n r = E @ (los_meas_noisy[idx] - los_pred)\n residuals.append(r)\n\n # IMPORTANT: Jacobian must be evaluated on un-normalized LOS vector\n U = unit_vector_jacobian(los_vec)\n # A is our state transition matrix for this time,\n # which maps changes in initial state to changes in state at this time.\n # We only care about the position part of the state\n # since the measurement is only a function of position.\n A = phi[idx][0:3, :]\n Hk = -E @ U @ A\n J_list.append(Hk)\n\n # OK, now do the least squares solve to get state update\n r = np.hstack(residuals)\n H = np.vstack(J_list)\n # Solve for state update using least squares\n dx = np.linalg.lstsq(H, r, rcond=None)[0]\n # Adjust the initial state estimate by the state update\n state0_est -= dx\n\nprint(\n f\"Initial State Error: {np.linalg.norm(state0_est_prior - state0_truth) / 1e3:.3f} km, {np.linalg.norm(state0_est_prior[3:6] - state0_truth[3:6]):.3f} m/s\"\n)\nprint(\n f\"Final State Error: {np.linalg.norm(state0_est - state0_truth) / 1e3:.3f} km, {np.linalg.norm(state0_est[3:6] - state0_truth[3:6]):.3f} m/s\"\n)"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}